	/*This program computes logit demand and simulates corresponding price changes assuming Bertrand competition with product differentiation*/

	#delimit;
	set matsize 10000;
	clear;
	clear matrix;
	mata: mata clear;
	local nbrands=7;
	set seed 070910432;

	*quietly{;
		set mem 250m;
		set more off;
		do "c:\data\restat_12_9\programs\csolve_do\csolve";
		do "c:\data\restat_12_9\programs\csolve_do\callf";
		do "c:\data\restat_12_9\programs\csolve_do\callf_setup";
		do "c:\data\restat_12_9\programs\csolve_do\antilogit_7_8";
		do "c:\data\restat_12_9\programs\csolve_do\anti_costs";
		do "c:\data\restat_12_9\programs\csolve_do\mcosts_antisyr";
		
		/*get top level elasticity*/
		mata: mata matuse "c:\data\restat_12_9/analysis_data/oil_top_level";
		
		mata: st_matrix("ols_top_b",ols_b_top[1,1]);
		mata: st_matrix("iv_top_b",iv_b_top[1,1]);
		
		mata: st_matrix("ols_top_mean",ols_b_top);
		mata: st_matrix("iv_top_mean",iv_b_top);
		mata: st_matrix("ols_top_vc",ols_vc_top);
		mata: st_matrix("iv_top_vc",iv_vc_top);

		
		
		scalar ols_top_b=ols_top_b[1,1];
		mata: ols_top_b=ols_b_top[1,1];
		
		scalar iv_top_b=iv_top_b[1,1];
		mata: iv_top_b=iv_b_top[1,1];
		
		use "c:\data\restat_12_9\analysis_data\hl_oillogitdata";
		qui tab BRAND2, gen(bra);
		qui tab REGION, gen(reg);
		qui tab date, gen(t);
		qui tab MONTH, gen(m);
		local reglist "reg2";
		forval x=3/10{;
			local reglist "`reglist' reg`x'";
		};
		
		su AVOL;
		scalar vt=r(N)*r(mean);
		forval x=1/`nbrands'{;
			su AVOL if bra`x'==1;
			scalar vshare`x'=r(N)*r(mean)/vt;
		};
		
		forval x=1/`nbrands'{;
			su price if bra`x'==1;
			scalar price`x'=r(mean);
		};
		scalar wavgp=vshare1*price1+vshare2*price2+vshare3*price3+vshare4*price4+vshare5*price5+vshare6*price6+vshare7*price7;
		forval x=1/`nbrands'{;
			forval y=2/10{;
				gen b`x'r`y'=reg`y'*bra`x';
			};
		};
		
		forval x=1/`nbrands'{;
			forval y=2/12{;
				gen m`x'r`y'=bra`x'*m`y';
			};
		};
		
		forval x=1/`nbrands'{;
			gen trend`x'=bra`x'*date;
		};
		
		
		forval x=1/4{;
			forval y=2/10{;
				local brandxreg "`brandxreg' b`x'r`y'";
			};
		};
		
		forval x=6/7{;
			forval y=2/10{;
				local brandxreg "`brandxreg' b`x'r`y'";
			};
		};
		
		forval x=1/4{;
			forval y=2/12{;
				local brandxmonth "`brandxmonth' m`x'r`y'";
			};
		};
		
		forval x=6/7{;
			forval y=2/12{;
				local brandxmonth "`brandxmonth' m`x'r`y'";
			};
		};
		
		local trends "trend6 trend7";
		forval x=1/4{;
			local trends "`trends' trend`x'";
		};
		
		local logit "price_plprice bra1 bra2 bra3 bra4 bra6 bra7 `reglist' t2 t3 t4 t5 t6 t7 t8 t9 t10 t11 t12 t13 t14 t15 t16 t17 t18 t19 t20 t21 t22 t23";
		local logit2 "price_plprice bra1 bra2 bra3 bra4 bra6 bra7 `brandxreg' `brandxmonth' `trends'";
		sort BRAND REGION date;
		
		
		sort BRAND REGION date;
		by BRAND REGION: gen T=_n;
		sort T REGION BRAND;
		bysort T: gen id=_n;
		tsset id T;
		
		noisily xi: newey2 lshare_plshare `logit'  if BRAND2~="PRIVATE", lag(1) nocons;
		noisily xi: newey2 lshare_plshare `logit2'  if BRAND2~="PRIVATE", lag(1) nocons;
		
		matrix ols_logitB=e(b);
		matrix ols_logitVC=e(V);
		scalar alpha=_b[price_plprice];
		matrix olsb=_b[price_plprice];
		
		
		forval x=1/`nbrands'{;
			forval y=1/`nbrands'{;
				sca olse`x'b`y'=-price`x'*alpha*vshare`x'+ols_top_b*vshare`x'*price`x'/wavgp;
			};
		};
		
		forval y=1/`nbrands'{;
			sca olse`y'b`y'=price`y'*alpha*(1-vshare`y')+ols_top_b*vshare`y'*price`y'/wavgp;
		};
		
		mat olselas=J(7,7,0);
		forval x=1/`nbrands'{;
			forval y=1/`nbrands'{;
				mat olselas[`x',`y']=olse`x'b`y';
			};
		};
		mat li olselas;
		
		
		noisily xi: newey2 lshare_plshare bra1 bra2 bra3 bra4 bra6 bra7  `reglist' t2 t3 t4 t5 t6 t7 t8 t9 t10 t11 t12 t13 t14 t15 t16 t17 t18 t19 t20 t21 t22 t23 (price_plprice=iprice) if BRAND2~="PRIVATE", nocons lag(1);
		noisily xi: newey2 lshare_plshare bra1 bra2 bra3 bra4 bra6 bra7 `brandxreg' `brandxmonth' `trends'  (price_plprice=iprice) if BRAND2~="PRIVATE", nocons lag(1);
		matrix iv_logitB=e(b);
		matrix iv_logitVC=e(V);
		
		scalar ivalpha=_b[price_plprice];
		matrix ivb=_b[price_plprice];
		
		
		forval x=1/`nbrands'{;
			forval y=1/`nbrands'{;
				sca ive`x'b`y'=-ivalpha*vshare`x'*price`x'+iv_top_b*vshare`x'*price`x'/wavgp;
			};
		};
		
		forval y=1/`nbrands'{;
			sca ive`y'b`y'=ivalpha*price`y'*(1-vshare`y')+iv_top_b*vshare`y'*price`y'/wavgp;
		};
		mat ivelas=J(7,7,0);
		forval x=1/`nbrands'{;
			forval y=1/`nbrands'{;
				mat ivelas[`x',`y']=ive`x'b`y';
			};
		};
		mat li ivelas;
		
 
		#delimit;

		mat vsbar=J(7,10,0);
		mat pbar=J(7,10,0);
		forval r=1/10{;
			sum AVOL if reg`r'==1;
			sca vt`r'=r(N)*r(mean);
		};
	
	
	
		forval r=1/10{;
			forval b=1/7{;
				qui su AVOL if reg`r'==1&bra`b'==1;
				mat vsbar[`b',`r']=r(N)*r(mean)/vt`r';
				qui su price if reg`r'==1&bra`b'==1;
				mat pbar[`b',`r']=r(mean);
			};	
		};		

		
		
	#delimit cr;
		mata:
			mata matuse c:/data/restat_12_9/analysis_data/oilhlpost_prices, replace
			pbar=st_matrix("pbar")
			vsbar=st_matrix("vsbar")
			postdirect=J(7,10,0)
			postdirect[1,.]=pbar[1,.]*(1+0.0141)	
			postdirect[2,.]=pbar[2,.]*(1-0.2060)	
			postdirect[3,.]=pbar[3,.]*(1-.0236)	
			postdirect[4,.]=pbar[4,.]*(1+0.0396)	
			postdirect[5,.]=pbar[5,.]*(1-0.0480)	
			postdirect[6,.]=pbar[6,.]*(1+0.0549)	
			postdirect[7,.]=pbar[7,.]*(1-0.0448)	
			
			
			pre_m=J(7,7,1)
			
			
			
			post_m=pre_m
			
			olslogit=J(`nbrands',10,0)
			ivlogit=J(`nbrands',10,0)
			ols_solver_flag=J(10,1,0)
			iv_solver_flag=J(10,1,0)
			olsb=st_matrix("olsb")
			ivb=st_matrix("ivb")
			olsb=-olsb
			ivb=-ivb
			ols_e=-ols_top_b
			iv_e=-iv_top_b
		
			
			for(i=1;i<=cols(pbar);i++){
					pinit=pbar[.,i]
					a=csolve(&anti(),pinit,1e-6,400,vsbar[.,i],pbar[.,i],olsb,ols_e,pre_m,post_m)
					csol=a[1..7,.]
					olslogit[.,i]=csol
					ols_solver_flag[i,1]=a[8,1]
			}
		
		olslogit=olslogit'
		ans=(olslogit-pbar'):/pbar'
		/*take only markets for which we could compute a solution*/
		ols_sim=100*10/(colsum(J(10,1,1)-!rowsum(ans)))*mean(ans,1)
		
			for(i=1;i<=cols(pbar);i++){
					pinit=pbar[.,i]
					a=csolve(&anti(),pinit,1e-6,400,vsbar[.,i],pbar[.,i],ivb,iv_e,pre_m,post_m)
					sol=a[1..7,.]
					ivlogit[.,i]=sol
					iv_solver_flag[i,1]=a[8,1]
			}
		ivlogit=ivlogit'
		ans=(ivlogit-pbar'):/pbar'
		iv_sim=100*10/(colsum(J(10,1,1)-!rowsum(ans)))*mean(ans,1)
		
		
		/*Pre-merger marginal costs, OLS*/
		pre_costs_ols=J(10,7,0)
		for(i=1;i<=cols(pbar);i++){
			pre_costs_ols[i,.]=anti_costs(vsbar[.,i],pbar[.,i],olsb,ols_e,pre_m)'
		}
		
			/*necessary marginal cost changes for observed price increase and simulations to match*/
			post_costs_ols=J(10,7,0)
			for(i=1;i<=cols(pbar);i++){
					init=pre_costs_ols[i,.]'
					a=csolve(&mcosts_antisyr(),init,1e-6,400,vsbar[.,i],postdirect[.,i],olsb,ols_e,pre_m,post_m)
					post_costs_ols[i,.]=a[1..7,.]'				
			}
			mcchange_ols=(post_costs_ols-pre_costs_ols):/pre_costs_ols
			permcchange_ols=100*10/(colsum(J(10,1,1)-!rowsum(mcchange_ols)))*mean(mcchange_ols,1)
		
		/*Pre-merger marginal costs, IV*/
		pre_costs_iv=J(10,7,0)
		for(i=1;i<=cols(pbar);i++){
			pre_costs_iv[i,.]=anti_costs(vsbar[.,i],pbar[.,i],ivb,iv_e,pre_m)'
		}
		
			post_costs_iv=J(10,7,0)
			for(i=1;i<=cols(pbar);i++){
					init=pre_costs_iv[i,.]'
					a=csolve(&mcosts_antisyr(),init,1e-6,400,vsbar[.,i],postdirect[.,i],ivb,iv_e,pre_m,post_m)
					post_costs_iv[i,.]=a[1..7,.]'				
			}
			mcchange_iv=(post_costs_iv-pre_costs_iv):/pre_costs_iv
			permcchange_iv=100*10/(colsum(J(10,1,1)-!rowsum(mcchange_iv)))*mean(mcchange_iv,1)
	end		
	

	
	

/***Now that we have our point estimates of costs and approximate simulated price changes, we draw from asymptotic approx of distributions
of estimators and recalculate for each draw.  will look at empirical distribution of these calculated costs and changes to do inference.****/	
#delimit;
drop _all;
/*number of bootstrap iterations*/
local its=1000;
drawnorm `logit2', n(`its') mean(ols_logitB) cov(ols_logitVC);
mkmat _all;


drop _all;
drawnorm e_topbs b2-b23, n(`its') mean(ols_top_mean) cov(ols_top_vc);
mkmat e_topbs;
drop _all;




use "c:\data\restat_12_9\analysis_data\hl_oillogitdata";

tab REGION, gen(reg);
tab BRAND2, gen(bra);
forval x=1/10{;
    mat bpbar`x'=J(`nbrands',`its',0);  
	mat bv`x'=J(`nbrands',`its',0);  
    mat bvsbar`x'=J(`nbrands',`its',0);  
};


forval x=1/10{;
    forval w=1/`its'{;
        preserve;
            keep if reg`x'==1;
            bsample;
            forval y=1/`nbrands'{;
                qui sum price if bra`y'==1;
                matrix bpbar`x'[`y',`w']=_result(3);
				qui sum AVOL if bra`y'==1;
                matrix bv`x'[`y',`w']=_result(3);				
            };
        restore;
    };
};
#delimit cr;
mata:
    bp1=st_matrix("bpbar1")
    bp2=st_matrix("bpbar2")
    bp3=st_matrix("bpbar3")
    bp4=st_matrix("bpbar4")
    bp5=st_matrix("bpbar5")
    bp6=st_matrix("bpbar6")
    bp7=st_matrix("bpbar7")
    bp8=st_matrix("bpbar8")
    bp9=st_matrix("bpbar9")
    bp10=st_matrix("bpbar10")
    

    ppoint=J(1,10,NULL)
    ppoint[1]=&bp1
    ppoint[2]=&bp2
    ppoint[3]=&bp3
    ppoint[4]=&bp4
    ppoint[5]=&bp5
    ppoint[6]=&bp6
    ppoint[7]=&bp7
    ppoint[8]=&bp8
    ppoint[9]=&bp9
    ppoint[10]=&bp10
        
	bv1=st_matrix("bv1")
    bv2=st_matrix("bv2")
    bv3=st_matrix("bv3")
    bv4=st_matrix("bv4")
    bv5=st_matrix("bv5")
    bv6=st_matrix("bv6")
    bv7=st_matrix("bv7")
    bv8=st_matrix("bv8")
    bv9=st_matrix("bv9")
    bv10=st_matrix("bv10") 
    
    vspoint=J(1,10,NULL)
    vspoint[1]=&bv1
    vspoint[2]=&bv2
    vspoint[3]=&bv3
    vspoint[4]=&bv4
    vspoint[5]=&bv5
    vspoint[6]=&bv6
    vspoint[7]=&bv7
    vspoint[8]=&bv8
    vspoint[9]=&bv9
    vspoint[10]=&bv10
     

    b=st_matrix("price_plprice")
	e=st_matrix("e_topbs")
end

    
	
set more off	
/*OLS Elasticities Here*/
local nbrands=7
mata:
    e1j=J(1000,7,0)
    e2j=J(1000,7,0)
    e3j=J(1000,7,0)
    e4j=J(1000,7,0)
    e5j=J(1000,7,0)
    e6j=J(1000,7,0)
    e7j=J(1000,7,0)
    
	

for(z=1;z<=1000;z++){
        /*makes price coefficient matrix for each draw*/
	bz=b[z,1]
	ez=e[z,1]
	
	/*need to get vs, p*/
	vx=J(7,10,0)
	px=J(7,10,0)

	for(j=1;j<=10;j++){           
		px[.,j]=(*ppoint[j])[.,z]
		vx[.,j]=(*vspoint[j])[.,z]
	}	            

	v=(rowsum(vx):/colsum(rowsum(vx)))'
	p=mean(px',1)

	
	elas=J(7,7,0)
	wavgp=v*p'
	for(m=1;m<=7;m++){
		for(n=1;n<=7;n++){
			elas[m,n]=-bz*v[m]*p[m]+ez*v[m]*p[m]/wavgp
		}
	}

	for(m=1;m<=7;m++){
		elas[m,m]=bz*p[m]*(1-v[m])+ez*v[m]*p[m]/wavgp
	}
		


e1j[z,1]=elas[1,1]
e1j[z,2]=elas[1,2]
e1j[z,3]=elas[1,3]
e1j[z,4]=elas[1,4]
e1j[z,5]=elas[1,5]
e1j[z,6]=elas[1,6]
e1j[z,7]=elas[1,7]


e2j[z,1]=elas[2,1]
e2j[z,2]=elas[2,2]
e2j[z,3]=elas[2,3]
e2j[z,4]=elas[2,4]
e2j[z,5]=elas[2,5]
e2j[z,6]=elas[2,6]
e2j[z,7]=elas[2,7]

e3j[z,1]=elas[3,1]
e3j[z,2]=elas[3,2]
e3j[z,3]=elas[3,3]
e3j[z,4]=elas[3,4]
e3j[z,5]=elas[3,5]
e3j[z,6]=elas[3,6]
e3j[z,7]=elas[3,7]

e4j[z,1]=elas[4,1]
e4j[z,2]=elas[4,2]
e4j[z,3]=elas[4,3]
e4j[z,4]=elas[4,4]
e4j[z,5]=elas[4,5]
e4j[z,6]=elas[4,6]
e4j[z,7]=elas[4,7]

e5j[z,1]=elas[5,1]
e5j[z,2]=elas[5,2]
e5j[z,3]=elas[5,3]
e5j[z,4]=elas[5,4]
e5j[z,5]=elas[5,5]
e5j[z,6]=elas[5,6]
e5j[z,7]=elas[5,7]

e6j[z,1]=elas[6,1]
e6j[z,2]=elas[6,2]
e6j[z,3]=elas[6,3]
e6j[z,4]=elas[6,4]
e6j[z,5]=elas[6,5]
e6j[z,6]=elas[6,6]
e6j[z,7]=elas[6,7]

e7j[z,1]=elas[7,1]
e7j[z,2]=elas[7,2]
e7j[z,3]=elas[7,3]
e7j[z,4]=elas[7,4]
e7j[z,5]=elas[7,5]
e7j[z,6]=elas[7,6]
e7j[z,7]=elas[7,7]

}

se_elas=J(7,7,0)

se_elas[1,.]=sqrt(mean(e1j:*e1j)-mean(e1j):*mean(e1j))				      
se_elas[2,.]=sqrt(mean(e2j:*e2j)-mean(e2j):*mean(e2j))				      
se_elas[3,.]=sqrt(mean(e3j:*e3j)-mean(e3j):*mean(e3j))
se_elas[4,.]=sqrt(mean(e4j:*e4j)-mean(e4j):*mean(e4j))
se_elas[5,.]=sqrt(mean(e5j:*e5j)-mean(e5j):*mean(e5j))				      
se_elas[6,.]=sqrt(mean(e6j:*e6j)-mean(e6j):*mean(e6j))				      
se_elas[7,.]=sqrt(mean(e7j:*e7j)-mean(e7j):*mean(e7j))				      

end	
		
	
	
	
	
	
	
	
	
	
	
	
	
	
mata:
	bslogit=J(1000,`nbrands',0)
	bslogitfocs=J(1000,`nbrands',0)
	flag=J(1000,10,0)	
					
	for(z=1;z<=1000;z++){
		bz=-b[z,1]
		ez=-e[z,1]		
        temp=J(10,7,0)

		for(j=1;j<=10;j++){
            	p=(*ppoint[j])[.,z]
	            pinit=(*ppoint[j])[.,z]
				v=(*vspoint[j])[.,z]
				vsbar=v:/colsum(v)
	            a=csolve(&anti(),pinit,1e-6,400,vsbar,p,bz,ez,pre_m,post_m)
      	        flag[z,j]=a[8,1]
				a=a[1..7,.]
            	
				a=(a-p):/p
	            temp[j,.]=a'
	      }
	bslogit[z,.]=100*10/(colsum(J(10,1,1)-!rowsum(temp)))*mean(temp,1)
	focs=J(10,7,0)
		for(j=1;j<=10;j++){
		    	p=(*ppoint[j])[.,z]	            
	            v=(*vspoint[j])[.,z]
				vsbar=v:/colsum(v)
				focs[j,.]=anti(postdirect[.,j],vsbar,p,bz,ez,pre_m,post_m)'
		}	
	bslogitfocs[z,.]=mean(focs)
	}
end
mata: st_matrix("bslogit", bslogit)
mata: st_matrix("bslogitfocs", bslogitfocs)
drop _all
svmat bslogit
svmat bslogitfocs
save "c:\data\restat_12_9\bootstrap_results\ols_oillogit_coll", replace





#delimit;
drop _all;
/*number of bootstrap iterations*/
local its=1000;
	local reglist "reg2";
	forval x=3/10{;
		local reglist "`reglist' reg`x'";
	};
local logit "price_plprice bra1 bra2 bra3 bra4 bra6 bra7 `brandxreg' `brandxmonth' `trends'";
drawnorm `logit', n(`its') mean(iv_logitB) cov(iv_logitVC);
mkmat _all;
drop _all;
drawnorm e_topbs b2-b23, n(`its') mean(iv_top_mean) cov(iv_top_vc);
mkmat e_topbs;
drop _all;

#delimit cr

/*IV Elasticities Here*/
mata:
	b=st_matrix("price_plprice")
	e=st_matrix("e_topbs")
    e1j=J(1000,7,0)
    e2j=J(1000,7,0)
    e3j=J(1000,7,0)
    e4j=J(1000,7,0)
    e5j=J(1000,7,0)
    e6j=J(1000,7,0)
    e7j=J(1000,7,0)
    
	

for(z=1;z<=1000;z++){
        /*makes price coefficient matrix for each draw*/
	bz=b[z,1]
	ez=e[z,1]
	
	/*need to get vs, p*/
	vx=J(7,10,0)
	px=J(7,10,0)

	for(j=1;j<=10;j++){           
		px[.,j]=(*ppoint[j])[.,z]
		vx[.,j]=(*vspoint[j])[.,z]
	}	            

	v=(rowsum(vx):/colsum(rowsum(vx)))'
	p=mean(px',1)

	
	elas=J(7,7,0)
	wavgp=v*p'
	for(m=1;m<=7;m++){
		for(n=1;n<=7;n++){
			elas[m,n]=-bz*v[m]*p[m]+ez*v[m]*p[m]/wavgp
		}
	}

	for(m=1;m<=7;m++){
		elas[m,m]=bz*p[m]*(1-v[m])+ez*v[m]*p[m]/wavgp
	}
		



e1j[z,1]=elas[1,1]
e1j[z,2]=elas[1,2]
e1j[z,3]=elas[1,3]
e1j[z,4]=elas[1,4]
e1j[z,5]=elas[1,5]
e1j[z,6]=elas[1,6]
e1j[z,7]=elas[1,7]


e2j[z,1]=elas[2,1]
e2j[z,2]=elas[2,2]
e2j[z,3]=elas[2,3]
e2j[z,4]=elas[2,4]
e2j[z,5]=elas[2,5]
e2j[z,6]=elas[2,6]
e2j[z,7]=elas[2,7]

e3j[z,1]=elas[3,1]
e3j[z,2]=elas[3,2]
e3j[z,3]=elas[3,3]
e3j[z,4]=elas[3,4]
e3j[z,5]=elas[3,5]
e3j[z,6]=elas[3,6]
e3j[z,7]=elas[3,7]

e4j[z,1]=elas[4,1]
e4j[z,2]=elas[4,2]
e4j[z,3]=elas[4,3]
e4j[z,4]=elas[4,4]
e4j[z,5]=elas[4,5]
e4j[z,6]=elas[4,6]
e4j[z,7]=elas[4,7]

e5j[z,1]=elas[5,1]
e5j[z,2]=elas[5,2]
e5j[z,3]=elas[5,3]
e5j[z,4]=elas[5,4]
e5j[z,5]=elas[5,5]
e5j[z,6]=elas[5,6]
e5j[z,7]=elas[5,7]

e6j[z,1]=elas[6,1]
e6j[z,2]=elas[6,2]
e6j[z,3]=elas[6,3]
e6j[z,4]=elas[6,4]
e6j[z,5]=elas[6,5]
e6j[z,6]=elas[6,6]
e6j[z,7]=elas[6,7]

e7j[z,1]=elas[7,1]
e7j[z,2]=elas[7,2]
e7j[z,3]=elas[7,3]
e7j[z,4]=elas[7,4]
e7j[z,5]=elas[7,5]
e7j[z,6]=elas[7,6]
e7j[z,7]=elas[7,7]

}

se_elas=J(7,7,0)

se_elas[1,.]=sqrt(mean(e1j:*e1j)-mean(e1j):*mean(e1j))				      
se_elas[2,.]=sqrt(mean(e2j:*e2j)-mean(e2j):*mean(e2j))				      
se_elas[3,.]=sqrt(mean(e3j:*e3j)-mean(e3j):*mean(e3j))
se_elas[4,.]=sqrt(mean(e4j:*e4j)-mean(e4j):*mean(e4j))
se_elas[5,.]=sqrt(mean(e5j:*e5j)-mean(e5j):*mean(e5j))				      
se_elas[6,.]=sqrt(mean(e6j:*e6j)-mean(e6j):*mean(e6j))				      
se_elas[7,.]=sqrt(mean(e7j:*e7j)-mean(e7j):*mean(e7j))				      

end	



#delimit cr;
mata:
    b=st_matrix("price_plprice")
	e=st_matrix("e_topbs")
    bslogit=J(1000,7,0)
	bslogitfocs=J(1000,7,0)

	for(z=1;z<=1000;z++){
		bz=-b[z,1]
		ez=-e[z,1]		
        temp=J(10,7,0)

		for(j=1;j<=10;j++){
            	p=(*ppoint[j])[.,z]
	            pinit=(*ppoint[j])[.,z]
				v=(*vspoint[j])[.,z]
				vsbar=v:/colsum(v)
	            a=csolve(&anti(),pinit,1e-6,400,vsbar,p,bz,ez,pre_m,post_m)
      	        a=a[1..7,.]
            	a=(a-p):/p
	            temp[j,.]=a'
	      }
	bslogit[z,.]=100*10/(colsum(J(10,1,1)-!rowsum(temp)))*mean(temp,1)
	
	focs=J(10,7,0)
		for(j=1;j<=10;j++){
		    	p=(*ppoint[j])[.,z]	            
	            v=(*vspoint[j])[.,z]
				vsbar=v:/colsum(v)
				focs[j,.]=anti(postdirect[.,j],vsbar,p,bz,ez,pre_m,post_m)'
		}	
	bslogitfocs[z,.]=mean(focs)
	}
end
mata: st_matrix("bslogit", bslogit)
mata: st_matrix("bslogitfocs", bslogitfocs)
drop _all
svmat bslogit
svmat bslogitfocs
save "c:\data\restat_12_9\bootstrap_results\iv_oillogit_coll", replace
